Influence of polarizability on metal oxide properties studied by molecular dynamics 

simulations 
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We have studied the dependence of metal oxide properties in molecular dynamics (MD) simula- 
tions on the polarizability of oxygen ions. We present studies of both liquid and crystalline structures 
of silica (SiC>2), magnesia (MgO) and alumina (AI2O3). For each of the three oxides, two separately 
optimized sets of force fields were used: (i) Long-range Coulomb interactions between oxide and 
metal ions combined with a short-range pair potential, (ii) Extension of force field (i) by adding 
polarizability to the oxygen ions. We show that while an effective potential of type (i) without 
polarizable oxygen ions can describe radial distributions and lattice constants reasonably well, po- 
tentials of type (ii) are required to obtain correct values for bond angles and the equation of state. 
The importance of polarizability for metal oxide properties decreases with increasing temperature. 
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I. INTRODUCTION 

Metal oxides are abundant in many technological ap- 
plications. Their excellent insulating, thermal isolating 
and heat resisting properties make them important com- 
ponents in microelectronics and semiconductor engineer- 
ing. They also play a crucial role in nanotcchnology 
and nanoscicncc. Modelling these systems in classical 
atomistic simulations requires effective potentials that 
describe the interaction between the ions with reason- 
able accuracy. There exists a wide selection of potential 
models with various degrees of sophistication, from sim- 
ple pair potentials to more intricate force fields. All of 
these have to cope with the long-ranged nature of the 
Coulomb interaction, which requires special treatment in 
large-scale simulations. In this work, we will show that 
certain properties of metal oxides cannot be adequately 
described in simulation with pure pair potentials. 

For crystalline silica (SiC^), Herzbach and Binder [1] 
already showed, that the polarizable force field proposed 
by Tangney and Scandolo [2] is superior to both the 
BKS [3] and the fluctuating-chargc DCG [4] pair poten- 
tial models. Recently, polarization effects on different 
properties of various molten fluorides, chlorides and ionic 
oxides were described by Salanne and Madden [5]. The 
authors illustrated the impact of polarizability in several 
cxemplarily selected systems of these material classes and 
predicted the general importance for structural and dy- 
namic properties. 

In this paper we present a systematic comparison of 
two types of force fields for three different oxides, which 
only differ by the presence of a polarizable term; the 
non-electrostatic and Coulomb terms have in each case 
the same functional form for both force fields. Hence, 
differences in molecular dynamics (MD) simulation ap- 
plications can be attributed to polarizability. To ob- 
tain a conclusion with a high degree of universality, we 
present MD simulation studies of both liquid and crys- 
talline structures of silica (SiC^), magnesia (MgO) and 



alumina (AI2O3). Simulations were performed with the 
MD code IMD [6, 7]. 

For each of the three metal oxides, we used two sepa- 
rately optimized force fields: 

(i) Short range interactions of Morse-Stretch (MS) 
form combined with Coulomb interactions between 
charged particles. 

(ii) Extension of (i) by adding polarizability to the oxy- 
gen ions according to the Tangney- Scandolo (TS) 
potential model [2]. 

The interactions between charges and/or induced 
dipoles are long ranged. To handle these electrostatic 
forces correctly and efficiently, we applied the Wolf sum- 
mation method [8] in the potential generation as well as 
in simulations. In several recent studies for metal oxides 
[9-11], linear scaling of the computational effort in the 
number of particles could be achieved by using the Wolf 
summation without significant loss of accuracy. 

Details of both TS model and Wolf summation are 
shown in Sec. II, where also the two force field types are 
described in detail. In Sec. Ill, whe present a system- 
atic comparison of the two types of force fields for silica 
(Si02), magnesia (MgO) and alumina (AI2O3). Discus- 
sion and conclusion are given in Sec. IV and Sec. V, re- 
spectively. 



II. FORCE FIELDS 

A. Generation 

All force fields were developed with the program potfit 
[12, 13], which generates effective interaction potentials 
solely from ab initio reference structures. The poten- 
tial parameters are optimized by matching the result- 
ing forces, energies, and stresses to corresponding first- 
principles values with the force matching method [14] . In 



2 





Dsi-Si 


Asi-o 


Do-o 


7Si-Si 


7Si-0 


70-0 




0.100 032 


0.100 250 


0.076 883 


11.009 449 


11.670 530 


7.505 632 




PSi-Si 


PSi-O 


po-o 


qsi 


go 






2.375 990 


2.073 780 


3.683 866 


1.799 475 


-0.899 738 






-DMg-Mg 


^Mg-O 


Do-o 


7Mg-Mg 


7m s -o 


7o-o 




0.038 258 


0.100 261 


0.065 940 


9.108 854 


10.405 694 


7.962 500 




PMg-Mg 


pMg-O 


po-o 


<7Mg 


go 






o.oo4 UUU 


0/117 QQQ 


UDU 




1 1 nn 79n 

-1.1UU 1 oU 






£>ai-ai 


£>ai-o 


Do-o 


7A1-A1 


7A1-0 


70-0 




0.002 164 


1.000 003 


0.000 018 


10.855 181 


7.617 923 


16.719 817 




Pai-ai 


Pai-o 


po-o 


<7ai 


qo 






5.517 666 


1.880 153 


6.609 171 


1.244 690 


-0.829 793 





TABLE I. Parameters for the potentials 



and <p A , given in IMD units set eV, A and amu. 



contrast to directly deriving charges and polarizabilities 
from first principles (as for example done in [15]) or fix- 
ing charges to their formal values (as for example done in 
[16]), these values are also to optimize in potfit in order 
to optimally reproduce forces, energies and stresses. It 
was already shown in [2, 10, 11], that such optimal force- 
fits yield highly accurate interaction potentials. This also 
implies that charges and polarizabilities obtained in this 
way are empirical quantities which need not correspond 
perfectly to a physical charge or polarizability. 

All reference data used in this study were obtained 
with the plane wave code VASP [17, 18]. For the opti- 
mization, a target function 



Z = w c Z c + w s Z s + Z{ 



(1) 



is minimized. Here, Z c , Z s and Z{ arc the contributions of 
the quadratic deviations of energies (e) , stresses (s) and 
forces (f), respectively. w e and w s are certain weights to 
balance the different amount of available data for each 
quantity. Root mean quare (RMS) errors, AFi (with I — 
c, s, f), are defined as proportional to the square root of 
the corresponding Zi . Their magnitudes are independent 
of weighting factors, number, and sizes of reference struc- 
tures. For the minimization of the potential parameters, 
a combination of a stochastic simulated annealing algo- 
rithm [19] and a conjugatc-gradicnt-likc deterministic al- 
gorithm [20] is used. Details of the whole optimization 
approach can be found in Rcf. 10 or 12. 



B. Wolf summation 

The electrostatic energies of a condensed system are 
commonly computed with the Ewald method [21], where 
the total Coulomb energy of a set of N ions, U qq , is de- 



composed into two terms U qq 
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r is the spacial coordinate. The short-ranged erfc term is 
summed up directly, while the smooth erf term is Fourier 
transformed and evaluated in reciprocal space. This re- 
stricts the technique to periodic systems. However, the 
main disadvantage is the scaling of the computational ef- 
fort with the number of particles in the simulation box, 
which increases as 0(N 3 / 2 ) [22], even for the optimal 
choice of the splitting parameter k. 

Wolf et al. [8] designed a method with linear scal- 
ing properties O(N) for Coulomb interactions. By tak- 
ing into account the physical properties of the system, 
the reciprocal-space term U qq is disregarded. In addi- 
tion, a continuous and smooth cutoff of the remaining 
screened Coulomb potential U qq (r.ij) = Qiqj erfc («rjj)r, "J 1 
is adopted at r c by shifting the potential so that it goes to 
zero smoothly in the first two derivatives at r = r c . The 
Wolf method for charges and its extension [9] to dipolar 
interactions is applied to all force fields the present article 
deals with. A detailed description of the Wolf summation 
with focus on its application to dipole contributions can 
be found in Refs. 9, 10, and 11, where it was successfully 
applied to silica, magnesia, and alumina. 



C. Type (i): MS + charges 



The potential (with fj, = S, M, A for silica, magnesia 
and alumina respectively) consists of a short range part 
of MS form, and a Coulomb interaction between charged 
particles. The MS interaction between two atoms i and 
j at positions r 2 ; and r, has the form 



Omk.o — Dij 



with 



exp[7y(l 



^)] 
Pij 



2exp[2^(l-^)] 
2 Pij 



(3) 



= I I , rij = Tj — Ti and the model pa- 
rameters D^, jij and py, which have to be optimized. 
The Coulomb interaction between two atoms is U qq ^j = 
Qi ( lj' r ij 1 - The charge q v of each atom type (with v = 
O, Si, Mg, Al for oxygen, silicon, magnesium and alu- 
minum respectively) is determined during the potential 
optimization under the constraint of charge neutrality. 
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The long-range Coulomb interactions are treated with 
the Wolf summation method both in force field gener- 
ation and simulation to obtain 4> m ,ij- The total inter- 
action is obtained by summing over all pairs of atoms: 

= <f>MS + ^qq = Yl(<f>MS,ij + 0qq,tf)- (4) 

i,j 

i>j 

With potfit, we generated effective interaction potentials 
</>° for silica, magnesia and alumina. Table I shows the 
corresponding parameters. 

D. Tangney-Scandolo dipole model 

In the TS model [2] , the polarizability a of the oxygen 
atoms is taken into account. The dipole moments depend 
on the local electric field of the surrounding charges and 
dipoles. Hence, a self-consistent iterative solution has to 
be found. In the TS approach, a dipole moment p" at 
position r*i in iteration step n consists of an induced part 
due to an electric field E{ri) and a short-range part pf R 
due to the short-range interactions between charges qi 
and qj . Following Rowley et at [23] , this contribution is 
given by 

p in = a J2 q Jpi fij ( nj ) (5) 

with 

fij{rij) = Cij J2^^e- b ^. (6) 

fe=0 

fij(rij) was introduced ad hoc to account for multipolc 
effects of nearest neighbors and is a function of very short 
range. b{j is the reciprocal of the length scale over which 
the short-range interaction comes into play, aj deter- 
mines amplitude and sign of this contribution to the in- 
duced moment, a, bij, and Cjj have to be optimized ad- 
ditionally. Together with the induced part, one obtains 

p? = ai E( ri ; {//' : }, ,.v. {r ? -} i=llJV ) + pf\ (7) 

where E(ri) is the electric field at position T*j, which 
is determined by the dipole moments pj in the previous 
iteration step. Due to its excellent performance and scal- 
ing properties (see Fig. 1 and Ref. 24), the TS model was 
used for more than 50 publications in the past ten years. 

E. Type (ii): (i) + dipoles 

Beside charges, electric dipole moments on each oxy- 
gen ion are taken into account using the TS model. In 
addition to the interaction between charges (0 qq ), the 
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TABLE II. RMS errors from the optimization of the force 
fields (f)^ and 0° (with fi = S, M, A for silica, magnesia and 
alumina respectively). 

Wolf summation method is also applied [9] to the inter- 
action charge-dipole (</> qp ) and dipole-dipole (4> PP )- This 
yields the total interaction 

4> = 4>MS + 0qq + 0qp + 0pp- (8) 

In this work, we use the polarizablc force fields for silica 
[10], magnesia [10] and alumina [11] presented in earlier 
publications. These potentials were also generated with 
potfit. 

It must again be emphasized that while the potentials 
</>° and cj)^ share the same functional form for short- 
range and Coulomb interactions, their potential param- 
eters were optimized individually. Otherwise, a compar- 
ison would only lead to the trivial result: A potential 
where some contributions are omitted is no longer accu- 
rate. 



F. RMS errors and scaling properties 

Although the potentials </> M have three more parame- 
ters than the </>° (13 compared to 10), they do not de- 
scribe the reference data significantly better than the <jp . 
This is illustrated by the RMS errors, that are shown in 
Table II. Indeed, there is a trend favoring the polarizable 
potentials: seven of nine RMS error are better in the case 
of And in the case of silica and magnesia, the RMS 
errors are - on average - smaller for (j>^ than for <fp . For 
alumina, however, it is the other way round. Thus, the 
RMS errors are only first indicators of the quality of the 
generated force field, but they are not able to denote for 
the practicability of a potential model. 

All presented force fields are pure pair term poten- 
tials. Hence, simulations are less expensive compared to 
other approaches with many-body potentials like a three- 
body interaction approach. Apart from that, considering 
polarizabiliy takes simulations time, because a new self- 
consistent solution for all dipole moments has to be found 
in each MD time step. Taking dipoles into account slows 
down simulation by a constant factor (in the case of silica 
about 2.6). The number of steps in the self-consistency 
loop is independent of the system size (Fig. 1). Both <^>° 
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FIG. 1. Scaling of computational effort with system size for 
0s and (f>s- 4>s also scales linear with system size using Wolf 
summation, but is slower by a factor of about 2.6 compared 
to S . The scaling properties of simulations of magnesia and 
alumina are the same as those of silica: M is slower than 0° 
by a factor independent of system size. 
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TABLE III. Lattice constants and cohesive energy per AI2O3 
unit of a-alumina at zero Kelvin obtained with a and <\r A 
compared with ab initio results and literature data. 



and M yield linear scaling with the number of particles 
due to the Wolf summation. A comparison of the Wolf 
performance with two mesh-based methods can be found 
in Ref. 10. 



III. RESULTS 

It is known from Ref. 10 that force fields may also 
yield qualitative results beyond the range for which they 
were optimized, however such applications beyond the 
optimization range should be closely verified. In the fol- 
lowing, we focus the tests on the range for the force fields 
were trained, but we also show results outside this zone 
to demonstrate the transferability of the potentials. 



a (A) e (A) Si-0 (A) Si-O-Si 

0s 4~98 5A7 L67 139 s 

0s [10] 5.15 5.50 1.65 148.5° 

Theory 4.97 [28] 5.39 [28] 1.61 [29] 145° [29] 



TABLE IV. Lattice constants, Si-0 bond length and Si-O-Si 
angle of Q-quartz at 300 K compared with theoretical studies. 



A. Microstructural properties 

First, the influence of polarizability on microstructural 
properties is illustrated. The radial distribution func- 
tions for liquid silica (4896 atoms) at 3000-6000 K are, 
in each case, evaluated for 100 snapshots taken out of 
100 ps MD runs at the given temperature. The aver- 
aged curves are given in Fig. 2. The curves obtained 
with 0g are similar to the curves of 0s- The existing 
slight deviations decrease with increasing temperature, 
so the polarizability is more important for lower temper- 
atures. For magnesia (see Fig. 3), the radial distribution 
functions show comparable behaviour: small deviations 
between <jP M and 4>Mi that decrease with temperature. 
Apparently, the radial distribution in high temperature 
oxide melts does not require polarizable oxide ions. 

A stronger influence of polarizability is observed on 
bond angles. Fig. 4 (left) depicts the oxygen centered 
angle distribution at 3000-6000 K in liquid silica and 
magnesia, respectively. Both 0g and 0^ overestimate 
the region of lower angles and underestimate the re- 
gion of higher angles. Potentials with electrostatic dipolc 
moments yields the correct shift of the distributions to 
slightly higher angles. Again, deviations decrease with 
increasing temperature. Although 0a and A were op- 
timized for low-temperature crystalline structures, we 
probed the behavior for liquid alumina at 3000 K. Fig. 4 
(right) shows the oxygen centered angle distribution in 
alumina compared to a recent ab initio study [25]. Al- 
though the trend of shifting angles to higher values by 
allowing for polarizability is not reproduced in this case, 
polarizability yields a curve which is in better agreement 
to ab initio data. These results coincide with Ref. 5, 
where the authors stated that polarization effects in ionic 
systems play an important role in determining bond an- 
gles. 

The electrostatic dipole moments also influence crys- 
talline structure parameters. The lattice constants of a- 
alumina are given in Table III. a yields an accurate 
agreement both with ab initio and experimental data, 
whereas with A the lattice constants are overestimated 
(in each case around 2% deviation). However, both sta- 
bilize the trigonal crystal structure. We also determined 
lattice constants, Si-O-Si angle and Si-0 bond length 
for a-quartz at 300 K (see Table IV), which is outside 
the optimization range for the 0s and 0g potentials. The 
average relative deviation of all parameters is for both 
potentials very similar (2.6% for 0g and 2.3% for 0g). 
On closer inspection, 0g yields more accurate lattice con- 
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FIG. 2. Normalized radial distribution functions for Si-Si, Si-O and O-O in liquid silica at 3000-6000 K. The solid (dashed) 
curves belong to 0s (0s). 
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FIG. 3. Normalized radial distribution functions for Mg-Mg, Mg-O and O-O in liquid magnesia at 3000-6000 K. The solid 
(dashed) curves belong to 0m (0m)- 



stants, whereas 0s better reproduces the Si-O-Si angle 
and Si-0 bond length. This is consistent with Ref. 5 and 
the results above concerning liquid metal oxides, where 
also the polarizability is more important for an improved 
description of bond angles rather than for atomic dis- 
tances. 



B. Thermodynamic properties 

For a-alumina, we also calculated the cohesive energy 
(see Table III). 0a coincides with ab initio and experi- 
mental results, A overestimates the cohesive energy (av- 
eraged deviation of 8.3%). This clear deviation shows 
that electrostatic dipole moments have to be taken into 
account, when probing macroscopical system properties. 

To investigate the influence of polarizability on other 
thermodynamic properties, we show the equation of state 



of liquid silica (3100 K, see Fig. 5) and magnesia (5000 K 
and 10 000 K, see Fig. 6) respectively. Pressures were ob- 
tained as averages along constant-volume MD runs of 10 
ps following 10 ps of equilibration. The curves obtained 
with 0s and 0m coincide with ab initio results as well 
as with experiment in the case of silica. The potentials 
0g and M , however, show a clear underestimation of 
the volume, which illustrates the need for polarizability. 
The insufficiency of M does not decrease with increasing 
temperature as in the case of microstructural properties. 
For the equation of state, polarizability has to be taken 
into account regardless of the simulation temperature. 



IV. DISCUSSION 

Apparently, the equation of state of liquid oxides shows 
the most significant difference between 0° and M ; the 
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FIG. 4. Left: Normalized oxygen centered angle distributions in liquid silica and magnesia at 3000-6000 K for 0s,M (solid 
curves) and 0s :M (dashed curves), yields a shoulder between around 80 and 130 degrees, whereas the band at around 
135-170 degrees is underestimated. Similarly, </>m yields a shoulder between around 65 and 95 degrees, whereas the band at 
around 110-150 degrees is underestimated. Deviations decrease with increasing temperature. Right: The oxygen centered 
angle distribution in liquid alumina is shown at 3000 K for 0a and <j>\ and compared to an ab initio study [25]. 
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FIG. 5. Equation of state of liquid silica at 3100 K for cj>s [10] 
and (f>s compared with experiment [30] and ab initio calcula- 
tions [31]. 
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FIG. 6. Equation of state of liquid magnesia at 5000 K (left 
three curves) and 10 000 K (right three curves) for <f>M and 
4>m compared with ab initio calculations [32]. 



potentials without polarizability seem to lack a signifi- 
cant contribution to the pressure. This also applies to 
the (non-polarizable) BKS force field [3], which underes- 
timates pressure at fixed specific volume by a comparable 
amount (cf. Ref. 2). The additional pressure in simula- 
tions with (f)^ can however not directly be attributed to 
dipolar interactions. An analysis of the virial showed 
that </> qp and r/> pp only contribute about 1.2% of the to- 
tal virial (and thus the pressure); the higher pressure for 



these systems results almost exclusively from stronger 
MS and Coulomb contributions to the virial. As the 
atomic forces arc described with comparable precision 
for both sets of potentials, this implies that the dipo- 
lar interaction is required to obtain correct pressures and 
forces simultaneously, especially as the polarizable force 
fields have higher absolute values of the atomic charges. 

When looking at the parameters of the force fields, it is 
noticable that in the 0°, the MS potentials are stronger 
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at smaller atomic distances: The absolute value of the 
MS strengh is higher and the stretch length p irj (cf. 
Eq. (3)) is shorter in the non-polarizable potentials. This 
seems to indicate that in this case, MS is required to 
describe the atomic interactions for nearest neighbours, 
while the dipolar interactions provide these contributions 
for the force fields with polarization. 

A better insight is uncovered by inspecting in detail, 
how accurate the reference ab initio forces and stresses 
are reproduced by the particular force fields: Fig. 7 de- 
picts scatter plots, where for each single quantity (force 
component, stress component) the value computed by 
the potential is plotted against its reference data value. 
Hence, for perfect matching a point is placed on the bi- 
secting line. As can be seen from the below graphs of 
Fig. 7 showing the force components of each single atom, 
both force field types (j)^ and <^>° yield distributions scat- 
tered around the bisecting line. The only difference be- 
tween ifin and (j)^ is how accurate the bisecting line is 
hit. In the top graphs, however, where the stress compo- 
nents of each configuration are depicted, clear deviations 
are uncovered: In the silica case, just a slightly worse 
matching of (j>° s compared to 4>s is visible. But <fP M pro- 
duces unnatural meanderings from the optimal matching 
at the left end of the graph (negative stress component 
values). The failure becomes even more apparent in the 
alumina case, where <jP A underestimates the stresses over 
the whole data set. To sum up, the scatter plots predict 
less accuracy for the <jp , when investigating system prop- 
erties with a strong dependence on stress and pressure. 

V. CONCLUSION 

In summary, we illustrated over a wide range the 
influence of polarizability on structural and thermody- 
namic properties in liquid and crystalline systems of silica 
(SiC>2), magnesia (MgO) and alumina (AI2O3), by com- 
paring two distinct potentials for each material. As the 
functional form (but not the parameters) of the short- 
range part and the Coulomb term are identical, the de- 



viations between the results of M over 0° could be asso- 
ciated with the additional dipole terms. 

We systematically investigated where the effects of 
electrostatic dipole moments are more important and 
how the impact arises from additional interaction mech- 
anisms. The strongest influence of polarizability is ob- 
served in macroscopic thermodynamic properties as the 
equation of state and the cohesive energy. Also on a 
microscopic scale, the role of dipoles is visible. However, 
the influences are more relevant for bond angle formation 
than for atomic distances in both liquid and crystalline 
structures. The influence always decreases with increas- 
ing temperature. 

Although the three presented metal oxides differ 
among each other in their stoichiometric configuration, 
the influence of polarizability is similar. Hence, it is ex- 
pected that conclusions can be extended to other metal 
oxide systems. The results are not estimated to be lim- 
ited to binary oxides as long as the same interaction 
mechanisms are dominant. At present, the polarizable 
force field approach is applied to yttrium doped zirco- 
nia [33], where a similar dependency of system proper- 
ties on polarizability is observed. Statements concerning 
systems with different interaction mechanisms such as 
hydrogen bonds in water may exceed the present study. 

Previous comparisons of different oxide potentials [1, 5] 
have already shown that polarizable force fields are su- 
perior in many aspects. In the present work, we demon- 
strate directly, for which system properties dipole inter- 
actions are important and where the influence is negligi- 
ble. Especially for the equation of state, simple pair po- 
tentials cannot reproduce the pressure-density relation- 
ship correctly. 
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